MLE fit for three component binding - simulated data

In this notebook we will see how well we can reproduce Kd of a non-fluorescent ligand from simulated experimental data with a maximum likelihood function.



In [1]:

    
import numpy as np
import matplotlib.pyplot as plt

from scipy import optimize
import seaborn as sns

%pylab inline









    



Populating the interactive namespace from numpy and matplotlib






    



:0: FutureWarning: IPython widgets are experimental and may change in the future.



In [2]:

    
#Competitive binding function
#This function and its assumptions are defined in greater detail in this notebook: 
## modelling-CompetitiveBinding-ThreeComponentBinding.ipynb
def three_component_competitive_binding(Ptot, Ltot, Kd_L, Atot, Kd_A):
    """
    Parameters
    ----------
    Ptot : float
        Total protein concentration
    Ltot : float
        Total tracer(fluorescent) ligand concentration
    Kd_L : float
        Dissociation constant of the fluorescent ligand
    Atot : float
        Total competitive ligand concentration
    Kd_A : float
        Dissociation constant of the competitive ligand
        
    Returns
    -------
    P : float
        Free protein concentration
    L : float
        Free ligand concentration
    A : float
        Free ligand concentration
    PL : float
        Complex concentration
    Kd_L_app : float
        Apparent dissociation constant of L in the presence of A
        
    Usage
    -----
    [P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot, Ltot, Kd_L, Atot, Kd_A)
    """
    Kd_L_app = Kd_L*(1+Atot/Kd_A)                                
    PL = 0.5 * ((Ptot + Ltot + Kd_L_app) - np.sqrt((Ptot + Ltot + Kd_L_app)**2 - 4*Ptot*Ltot))  # complex concentration (uM)
    P = Ptot - PL; # free protein concentration in sample cell after n injections (uM)                                                                                                                                                                                                                          
    L = Ltot - PL; # free tracer ligand concentration in sample cell after n injections (uM)
    A = Atot - PL; # free competitive ligand concentration in sample cell after n injections (uM)
    return [P, L, A, PL, Kd_L_app]



In [3]:

    
#Let's define our parameters
Kd = 3800e-9 # M
Kd_Competitor = 3000e-9 # M
Ptot = 1e-9 * np.ones([12],np.float64) # M
Ltot = 20.0e-6 / np.array([10**(float(i)/2.0) for i in range(12)]) # M
L_Competitor =  10e-6 # M



In [4]:

    
[P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot, Ltot, Kd, L_Competitor, Kd_Competitor)
#usin
[P_base, L_base, A_base, PL_base, Kd_L_app_base] = three_component_competitive_binding(Ptot, Ltot, Kd, 0, Kd_Competitor)



In [5]:

    
# y will be complex concentration
# x will be total ligand concentration
plt.title('Competition assay')
plt.semilogx(Ltot,PL_base,'green', marker='o', label = 'Fluorescent Ligand')
plt.semilogx(Ltot,PL,'cyan', marker='o', label = 'Fluorescent Ligand + Competitor')
plt.xlabel('$[L]_{tot}$')
plt.ylabel('$[PL]$')
plt.legend(loc=0);



In [6]:

    
#What if we change our Kd's a little
#Kd_Gef_Abl
Kd = 480e-9 # M
#Kd_Ima_Abl 
Kd_Competitor = 21.0e-9 # M



In [7]:

    
[P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot, Ltot, Kd, L_Competitor, Kd_Competitor)
#using _base as a subscript to define when we have no competitive ligand
[P_base, L_base, A_base, PL_base, Kd_L_app_base] = three_component_competitive_binding(Ptot, Ltot, Kd, 0, Kd_Competitor)



In [8]:

    
# y will be complex concentration
# x will be total ligand concentration
plt.title('Competition assay')
plt.semilogx(Ltot,PL_base,'green', marker='o', label = 'Fluorescent Ligand')
plt.semilogx(Ltot,PL,'cyan', marker='o', label = 'Fluorescent Ligand + Competitor')
plt.xlabel('$[L]_{tot}$')
plt.ylabel('$[PL]$')
plt.legend(loc=0);

Now make this a fluorescence experiment.



In [9]:

    
# Making max 400 relative fluorescence units, and scaling all of PL to that
npoints = len(Ltot)
sigma = 10.0 # size of noise
F_i = (400/1e-9)*PL + sigma * np.random.randn(npoints)
F_i_base = (400/1e-9)*PL_base + sigma * np.random.randn(npoints)



In [10]:

    
# y will be complex concentration
# x will be total ligand concentration
plt.title('Competition assay')
plt.semilogx(Ltot,F_i_base,'green', marker='o', label = 'Fluorescent Ligand')
plt.semilogx(Ltot,F_i,'cyan', marker='o', label = 'Fluorescent Ligand + Competitor')
plt.xlabel('$[L]_{tot}$')
plt.ylabel('$Fluorescence$')
plt.legend(loc=0);



In [11]:

    
#And makeup an F_L
F_L = 0.3

First let's see if we can find Kd of our fluorescent ligand from the three component binding model, if we know there's no competitive ligand



In [12]:

    
# This function fits Kd_L when L_Competitor is 0
def find_Kd_from_fluorescence_base(params):
    
    [F_background, F_PL, Kd_L] = params
    
    N = len(Ltot)
    Fmodel_i = np.zeros([N])
    
    for i in range(N):
        [P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot[0], Ltot[i], Kd_L, 0, Kd_Competitor)
        Fmodel_i[i] = (F_PL*PL + F_L*L) + F_background
    
    return Fmodel_i



In [ ]:



In [13]:

    
initial_guess = [1,400/1e-9,3800e-9]
prediction = find_Kd_from_fluorescence_base(initial_guess)

plt.semilogx(Ltot,prediction,color='k')
plt.semilogx(Ltot,F_i_base, 'o')
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend();









    



/Users/hansons/anaconda/lib/python2.7/site-packages/matplotlib/axes/_axes.py:519: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.
  warnings.warn("No labelled objects found. "



In [14]:

    
def sumofsquares(params):
    prediction = find_Kd_from_fluorescence_base(params)
    return np.sum((prediction - F_i_base)**2)



In [15]:

    
initial_guess = [0,3E11,2000E-9]
fit = optimize.minimize(sumofsquares,initial_guess,method='Nelder-Mead')
print "The predicted parameters are", fit.x









    



The predicted parameters are [ -2.74841327e+00   4.03350513e+11   4.45151221e-07]



In [16]:

    
fit_prediction = find_Kd_from_fluorescence_base(fit.x)



In [17]:

    
plt.semilogx(Ltot,fit_prediction,color='k')
plt.semilogx(Ltot,F_i_base, 'o')
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend();



In [18]:

    
Kd_L_MLE = fit.x[2]



In [19]:

    
def Kd_format(Kd):
    if (Kd < 1e-12):
        Kd_summary = "Kd = %.1f nM " % (Kd/1e-15)
    elif (Kd < 1e-9):
        Kd_summary = "Kd = %.1f pM " % (Kd/1e-12)
    elif (Kd < 1e-6):
        Kd_summary = "Kd = %.1f nM " % (Kd/1e-9)
    elif (Kd < 1e-3):
        Kd_summary = "Kd = %.1f uM " % (Kd/1e-6)
    elif (Kd < 1):
        Kd_summary = "Kd = %.1f mM " % (Kd/1e-3)
    else:
        Kd_summary = "Kd = %.3e M " % (Kd)
    
    return Kd_summary



In [20]:

    
Kd_format(Kd_L_MLE)









    Out[20]:





'Kd = 445.2 nM '



In [21]:

    
delG_summary = "delG = %s kT" %np.log(Kd_L_MLE)
delG_summary









    Out[21]:





'delG = -14.6248517897 kT'

Okay, cool now let's try to fit our competitive ligand



In [22]:

    
# This function fits Kd_A when Kd_L already has an estimate

def find_Kd_from_fluorescence_competitor(params):
    
    [F_background, F_PL, Kd_Competitor] = params
    
    N = len(Ltot)
    Fmodel_i = np.zeros([N])
    
    for i in range(N):
        [P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot[0], Ltot[i], Kd_L_MLE, L_Competitor, Kd_Competitor)
        Fmodel_i[i] = (F_PL*PL + F_L*L) + F_background
    
    return Fmodel_i



In [23]:

    
initial_guess = [0,400/1e-9,3800e-9]
prediction = find_Kd_from_fluorescence_competitor(initial_guess)

plt.semilogx(Ltot,prediction,color='k')
plt.semilogx(Ltot,F_i, 'o')
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend();



In [24]:

    
def sumofsquares(params):
    prediction = find_Kd_from_fluorescence_competitor(params)
    return np.sum((prediction - F_i)**2)



In [25]:

    
initial_guess = [0,3E11,2000E-9]
fit_comp = optimize.minimize(sumofsquares,initial_guess,method='Nelder-Mead')
print "The predicted parameters are", fit_comp.x









    



The predicted parameters are [ -6.93274088e-04   2.12842582e+11   6.40589610e-08]



In [26]:

    
fit_prediction_comp = find_Kd_from_fluorescence_competitor(fit_comp.x)



In [27]:

    
plt.semilogx(Ltot,fit_prediction_comp,color='k')
plt.semilogx(Ltot,F_i, 'o')
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend();



In [28]:

    
Kd_A_MLE = fit_comp.x[2]



In [29]:

    
Kd_format(Kd_A_MLE)









    Out[29]:





'Kd = 64.1 nM '



In [30]:

    
delG_summary = "delG = %s kT" %np.log(Kd_A_MLE)
delG_summary









    Out[30]:





'delG = -16.5634619124 kT'

Let's plot all these results with each other and see how we did



In [31]:

    
plt.semilogx(Ltot,fit_prediction_comp,color='k')
plt.semilogx(Ltot,F_i, 'o')
plt.axvline(Kd_A_MLE,color='k',label='%s (MLE)'%Kd_format(Kd_A_MLE))
plt.axvline(Kd_Competitor,color='b',label='%s'%Kd_format(Kd_Competitor))

plt.semilogx(Ltot,fit_prediction,color='k')
plt.semilogx(Ltot,F_i_base, 'o')
plt.axvline(Kd_L_MLE,color='k',label='%s (MLE)'%Kd_format(Kd_L_MLE))
plt.axvline(Kd,color='g',label='%s'%Kd_format(Kd))
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend(loc=0);



In [32]:

    
#Awesome



In [ ]: